%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This file is part of the book
%%
%% Algorithmic Graph Theory
%% http://code.google.com/p/graph-theory-algorithms-book/
%%
%% Copyright (C) 2010 Nathann Cohen <nathann.cohen@gmail.com>
%%
%% See the file COPYING for copying conditions.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\chapter{Graph Problems and Their LP Formulations}

This document is meant as an explanation of several graph theoretical functions defined in Sage's Graph Library (http://www.sagemath.org/), which use Linear Programming to solve optimization of existence problems.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Maximum average degree}
\label{lp:mad}

The average degree of a graph $G$ is defined as $ad(G) = \frac {2|E(G)|}{|V(G)|}$. The maximum average degree of $G$ is meant to represent its densest part, and is formally defined as : $$mad(G) = \max_{H\subseteq G}ad(H)$$

Even though such a formulation does not show it, this quantity can be computed in polynomial time through Linear Programming. Indeed, we can think of this as a simple flow problem defined on a bipartite graph. Let $D$ be a directed graph whose vertex set we first define as the disjoint union of $E(G)$ and $V(G)$. We add in $D$ an edge between $(e,v)\in E(G)\times V(G)$ if and only if $v$ is one of $e$'s endpoints. Each edge will then have a flow of 2 (through the addition in $D$ of a source and the necessary edges) to distribute among its two endpoints. We then write in our linear program the constraint that each vertex can absorb a flow of at most $z$ (add to $D$ the necessary sink and the edges with capacity $z$).

Clearly, if $H\subseteq G$ is the densest subgraph in $G$, its $|E(H)|$ edges will send a flow of $2|E(H)|$ to their $|V(H)|$ vertices, such a flow being feasible only if $z\geq \frac {2|E(H)|}{|V(H)|}$. An elementary application of the max-flow/min-cut theorem, or of Hall's bipartite matching theorem shows that such a value for $z$ is also sufficient. This LP can thus let us compute the Maximum Average Degree of the graph.

{\bf Sage method :}\verb! Graph.maximum_average_degree()!

{\bf LP Formulation :}
\begin{itemize}
\item Minimize : $z$\\
\item Such that :
  \begin{itemize}
  \item     a vertex can absorb at most $z$
    $$\forall v\in V(G), \sum_{e\in E(G)\atop e\sim v}x_{e,v}\leq z$$
\item each edge sends a flow of 2
  $$\forall e=uv\in E(G), x_{e,u} + x_{e,u} = 2$$
  \end{itemize}
  \item $x_{e,v}$ real positive variable
\end{itemize}

Here is the corresponding Sage code:

\begin{lstlisting}
sage: g = graphs.PetersenGraph()
sage: p = MixedIntegerLinearProgram(maximization = False)
sage: x = p.new_variable( dim = 2 )

sage: p.set_objective(p['z'])

sage: for v in g:
...     p.add_constraint( sum([ x[u][v] for u in g.neighbors(v) ]) <= p['z'] )

sage: for u,v in g.edges(labels = False):
...     p.add_constraint( x[u][v] + x[v][u] == 2 )

sage: p.solve()
3.0
\end{lstlisting}

{\bf REMARK : } In many if not all the other LP formulations, this Linear Program is used as a constraint. In those problems, we are always at some point looking for a subgraph $H$ of $G$ such that $H$ does not contain any cycle. The edges of $G$ are in this case variables, whose value can be equal to $0$ or $1$ depending on whether they belong to such a graph $H$. Based on the observation that the Maximum Average Degree of a tree on $n$ vertices is exactly its average degree ($=2-2/n<1$), and that any cycles in a graph ensures its average degree is larger than $2$, we can then set the constraint that $z\leq 2-\frac 2 {|V(G)|}$. This is a handy way to write in LP the constraint that ``the set of edges belonging to $H$ is acyclic''. For this to work, though, we need to ensure that the variables corresponding to our edges are binary variables.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Traveling Salesman Problem}
Given a graph $G$ whose edges are weighted by a function $w:E(G) \to \R$, a solution to the $TSP$ is a Hamiltonian (spanning) cycle whose weight (the sum of the weight of its edges) is minimal. It is easy to define both the objective and the constraint that each vertex must have exactly two neighbors, but this could produce solutions such that the set of edges define the disjoint union of several cycles. One way to formulate this linear program is hence to add the constraint that, given an arbitrary vertex $v$, the set $S$ of edges in the solution must contain no cycle in $G-v$, which amounts to checking that the set of edges in $S$ no adjacent to $v$ is of maximal average degree strictly less than 2, using the remark from section \ref{mad}.

We will then, in this case, define variables representing the edges included in the solution, along with variables representing the weight that each of these edges will send to their endpoints.

{\bf LP Formulation :}
\begin{itemize}
\item Minimize $$\sum_{e\in E(G)}w(e) b_e$$\\
\item Such that :
  \begin{itemize}
  \item Each vertex is of degree 2
    $$\forall v\in V(G), \sum_{e\in E(G)\atop e\sim v}b_e = 2$$
  \item No cycle disjoint from a special vertex $v^*$
    \begin{itemize}
    \item Each edge sends a flow of 2 if it is taken
      $$\forall e=uv\in E(G-v^*), x_{e,u} + x_{e,v} = 2b_e$$
    \item Vertices receive strictly less than 2
      $$\forall v\in V(G-v^*), \sum_{e\in E(G)\atop e\sim v}x_{e,v}\leq 2-\frac 2 {|V(G)|}$$
    \end{itemize}

  \end{itemize}
\item Variables
  \begin{itemize}
  \item $x_{e,v}$ real positive variable (flow sent by the edge)
  \item $b_e$ binary (is the edge in the solution ?)

  \end{itemize}
\end{itemize}

{\bf Sage method : }\verb!Graph.traveling_salesman_problem()!

Here is the corresponding Sage corresponding to a simpler case -- looking for an Hamiltonian cycle in a graph:

\begin{lstlisting}
sage: g = graphs.GridGraph([4,4])
sage: p = MixedIntegerLinearProgram(maximization = False)

sage: f = p.new_variable()
sage: r = p.new_variable()

sage: eps = 1/(2*Integer(g.order()))
sage: x = g.vertex_iterator().next()

sage: # reorders the edge as they can appear in the two different ways
sage: R = lambda x,y : (x,y) if x < y else (y,x)

sage: # All the vertices have degree 2
sage: for v in g:
...     p.add_constraint( sum([ f[R(u,v)] for u in g.neighbors(v)]) == 2)

sage: # r is greater than f
sage: for u,v in g.edges(labels = None):
...     p.add_constraint( r[(u,v)] + r[(v,u)] - f[R(u,v)] >= 0)

sage: # no cycle which does not contain x
sage: for v in g:
...     if v != x:
...          p.add_constraint( sum([ r[(u,v)] for u in g.neighbors(v)]) <= 1-eps)

sage: p.set_objective(None)
sage: p.set_binary(f)

sage: p.solve()                                    # optional - GLPK,CBC,CPLEX
0.0

sage: # We can now build the solution
sage: # found as a graph

sage: f = p.get_values(f)                          # optional - GLPK,CBC,CPLEX
sage: tsp = Graph()                                # optional - GLPK,CBC,CPLEX
sage: for e in g.edges(labels = False):            # optional - GLPK,CBC,CPLEX
...       if f[R(e[0],e[1])] == 1:                 # optional - GLPK,CBC,CPLEX
...          tsp.add_edge(e)                       # optional - GLPK,CBC,CPLEX

sage: tsp.is_regular(k=2) and tsp.is_connected()   # optional - GLPK,CBC,CPLEX
True
sage: tsp.order() == g.order()                     # optional - GLPK,CBC,CPLEX
True
\end{lstlisting}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Edge-disjoint spanning trees}

This problem is polynomial by a result from Edmonds. Obviously, nothing ensures the following formulation is a polynomial algorithm as it contains many integer variables, but it is still a short practical way to solve it.

This problem amounts to finding, given a graph $G$ and an integer $k$, edge-disjoint spanning trees $T_1, \dots, T_k$ which are subgraphs of $G$. In this case, we will chose to define a spanning tree as an acyclic set of $|V(G)|-1$ edges.

{\bf Sage method : }\verb!Graph.edge_disjoint_spanning_trees()!

{\bf LP Formulation :}
\begin{itemize}
\item Maximize : nothing
\item Such that :
  \begin{itemize}
  \item An edge belongs to at most one set
    $$\forall e\in E(G), \sum_{i\in [1,\dots,k]} b_{e,k} \leq 1$$
  \item Each set contains $|V(G)|-1$ edges
    $$\forall i\in [1,\dots,k], \sum_{e\in E(G)} b_{e,k} = |V(G)|-1$$
  \item No cycles
    \begin{itemize}
    \item In each set, each edge sends a flow of 2 if it is taken
      $$\forall i\in [1,\dots,k], \forall e=uv\in E(G), x_{e,k,u} + x_{e,k,u} = 2b_{e,k}$$
    \item Vertices receive strictly less than 2
      $$\forall i\in [1,\dots,k], \forall v\in V(G), \sum_{e\in E(G)\atop e\sim v}x_{e,k,v}\leq 2-\frac 2 {|V(G)|}$$
    \end{itemize}
  \end{itemize}
\item Variables
  \begin{itemize}
  \item $b_{e,k}$ binary (is edge $e$ in set $k$ ?)
  \item $ x_{e,k,u}$ positive real (flow sent by edge $e$ to vertex $u$ in set $k$)
  \end{itemize}
\end{itemize}

Here is the corresponding Sage code:

\begin{lstlisting}
sage: g = graphs.RandomGNP(40,.6)
sage: p = MixedIntegerLinearProgram()
sage: colors = range(2)

sage: # Sort an edge
sage: S = lambda (x,y) : (x,y) if x<y else (y,x)

sage: edges = p.new_variable(dim = 2)
sage: r_edges = p.new_variable(dim = 2)

sage: # An edge belongs to at most one tree
sage: for e in g.edges(labels=False):
...       p.add_constraint(sum([edges[j][S(e)] for j in colors]), max = 1)

sage: for j in colors:
...       # each color class has g.order()-1 edges
...       p.add_constraint(
...         sum([edges[j][S(e)] for e in g.edges(labels=None)])
...         >= g.order()-1)
...       # Each vertex is in the tree
...       for v in g.vertices():
...            p.add_constraint(
...              sum([edges[j][S(e)] for e in g.edges_incident(v, labels=None)])
...              >= 1)
...       # r_edges is larger than edges
...       for u,v in g.edges(labels=None):
...            p.add_constraint(
...               r_edges[j][(u,v)] + r_edges[j][(v, u)]
...               == edges[j][S((u,v))] )

sage: # no cycles
sage: epsilon =  (3*Integer(g.order()))**(-1)
sage: for j in colors:
...      for v in g:
...          p.add_constraint(
...             sum([r_edges[j][(u,v)] for u in g.neighbors(v)])
...             <= 1-epsilon)

sage: p.set_binary(edges)
sage: p.set_objective(None)
sage: p.solve()                                     # optional - GLPK,CBC,CPLEX
0.0


sage: # We can now build the solution
sage: # found as a list of trees

sage: edges = p.get_values(edges)                   # optional - GLPK,CBC,CPLEX
sage: trees = [Graph() for c in colors]             # optional - GLPK,CBC,CPLEX

sage: for e in g.edges(labels = False):             # optional - GLPK,CBC,CPLEX
...       for c in colors:                          # optional - GLPK,CBC,CPLEX
...          if round(edges[c][S(e)]) == 1:         # optional - GLPK,CBC,CPLEX
...              trees[c].add_edge(e)               # optional - GLPK,CBC,CPLEX

sage: all([ trees[j].is_tree() for j in colors ])   # optional - GLPK,CBC,CPLEX
True
\end{lstlisting}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Steiner tree}

See Trietsch~\cite{Trietsch1984} for a relationship between Steiner
trees and Euler's problem of polygon division.
Finding  a spanning tree  in a Graph $G$ can be done in linear time, whereas computing a Steiner Tree is NP-hard. The goal is in this case, given a graph, a weight function $w:E(G) \to \R$ and a set $S$ of vertices, to find the tree of minimum cost connecting them all together. Equivalently, we will be looking for an acyclic subgraph $H$of $G$ containing $|V(H)|$ vertices and $|E(H)|=|V(H)|-1$ edges, which contains each vertex from $S$

{\bf LP Formulation :}
\begin{itemize}
\item Minimize : $$\sum_{e\in E(G)}w(e) b_e$$\\
\item Such that :
  \begin{itemize}
  \item Each vertex from $S$ is in the tree
    $$\forall v \in S, \sum_{e\in E(G)\atop e\sim v}b_e \geq 1$$
  \item $c$ is equal to 1 when a vertex $v$ is in the tree
    $$\forall v \in V(G), \forall e\in E(G), e\sim v, b_e\leq c_v$$
  \item The tree contains $|V(H)|$ vertices and $|E(H)|=|V(H)|-1$ edges
    $$\sum_{v\in G}c_v - \sum_{e\in E(G)}b_e = 1$$
  \item No Cycles
    \begin{itemize}
    \item Each edge sends a flow of 2 if it is taken
      $$\forall e=uv\in E(G), x_{e,u} + x_{e,u} = 2b_{e,k}$$
    \item Vertices receive strictly less than 2
      $$\forall v\in V(G), \sum_{e\in E(G)\atop e\sim v}x_{e,v}\leq 2-\frac 2 {|V(G)|}$$
    \end{itemize}

  \end{itemize}
\item Variables :
  \begin{itemize}
  \item $b_e$  binary (is $e$ in the tree ?)
  \item $c_v$ binary (does the tree contain $v$ ?)
  \item $x_{e,v}$ real positive variable (flow sent by the edge)
  \end{itemize}
\end{itemize}


{\bf Sage method :}\verb! Graph.steiner_tree()!

Here is the corresponding Sage code:

\begin{lstlisting}
sage: g = graphs.GridGraph([10,10])
sage: vertices = [(0,2),(5,3)]

sage: from sage.numerical.mip import MixedIntegerLinearProgram
sage: p = MixedIntegerLinearProgram(maximization = False)

sage: # Reorder an edge
sage: R = lambda (x,y) : (x,y) if x<y else (y,x)

sage: # edges used in the Steiner Tree
sage: edges = p.new_variable()
sage: # relaxed edges to test for acyclicity
sage: r_edges = p.new_variable()
sage: # Whether a vertex is in the Steiner Tree
sage: vertex = p.new_variable()

sage: #  Which vertices are in the tree ?
sage: for v in g:
...         for e in g.edges_incident(v, labels=False):
...             p.add_constraint(vertex[v] - edges[R(e)], min = 0)

sage: # We must have the given vertices in our tree
sage: for v in vertices:
...         p.add_constraint(
...           sum([edges[R(e)] for e in g.edges_incident(v,labels=False)]
...            == 1)

sage: # The number of edges is equal to the number of vertices in our tree minus 1
sage: p.add_constraint(
...     sum([vertex[v] for v in g])
...     - sum([edges[R(e)] for e in g.edges(labels=None)])
...     == 1)

sage: # There are no cycles in our graph
sage: for u,v in g.edges(labels = False):
...         p.add_constraint(
...           r_edges[(u,v)]+ r_edges[(v,u)] - edges[R((u,v))]
...           <= 0 )

sage: eps = 1/(5*Integer(g.order()))

sage: for v in g:
...         p.add_constraint(sum([r_edges[(u,v)] for u in g.neighbors(v)]), max = 1-eps)

sage: p.set_objective(sum([edges[R(e)] for e in g.edges(labels = False)]))
sage: p.set_binary(edges)
sage: p.solve()                                    # optional - GLPK,CBC,CPLEX
6.0

sage: # We can now build the solution
sage: # found as a tree

sage: edges = p.get_values(edges)                  # optional - GLPK,CBC,CPLEX
sage: st = Graph()                                 # optional - GLPK,CBC,CPLEX
sage: st.add_edges(
...     [e for e in g.edges(labels = False)
...     if edges[R(e)] == 1])                      # optional - GLPK,CBC,CPLEX
sage: st.is_tree()                                 # optional - GLPK,CBC,CPLEX
True
sage: all([v in st for v in vertices])             # optional - GLPK,CBC,CPLEX
True
\end{lstlisting}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Linear arboricity}
The linear arboricity of a graph $G$ is the least number $k$ such that the edges of $G$ can be partitioned into $k$ classes, each of them being a forest of paths (the disjoints union of paths -- trees of maximal degree 2). The corresponding LP is very similar to the one giving edge-disjoint spanning trees

{\bf LP Formulation :}
\begin{itemize}
\item Maximize : nothing
\item Such that :
  \begin{itemize}
  \item An edge belongs to exactly one set
    $$\forall e\in E(G), \sum_{i\in [1,\dots,k]} b_{e,k}  = 1$$
  \item Each class has maximal degree 2
    $$\forall v\in V(G), \forall i\in [1,\dots,k], \sum_{e\in E(G)\atop e\sim v} b_{e,k}\leq 2$$
  \item No cycles
    \begin{itemize}
    \item In each set, each edge sends a flow of 2 if it is taken
      $$\forall i\in [1,\dots,k], \forall e=uv\in E(G), x_{e,k,u} + x_{e,k,v} = 2b_{e,k}$$
    \item Vertices receive strictly less than 2
      $$\forall i\in [1,\dots,k], \forall v\in V(G), \sum_{e\in E(G)\atop e\sim v}x_{e,k,v}\leq 2-\frac 2 {|V(G)|}$$
    \end{itemize}
  \end{itemize}
\item Variables
  \begin{itemize}
  \item $b_{e,k}$ binary (is edge $e$ in set $k$ ?)
  \item $ x_{e,k,u}$ positive real (flow sent by edge $e$ to vertex $u$ in set $k$)
  \end{itemize}
\end{itemize}

{\bf Sage method :}\verb! sage.graphs.graph_coloring.linear_arboricity()!

Here is the corresponding Sage code :

\begin{lstlisting}
sage: g = graphs.GridGraph([4,4])
sage: k = 2
sage: p = MixedIntegerLinearProgram()

sage: # c is a boolean value such that c[i][(u,v)] = 1
sage: # if and only if (u,v) is colored with i
sage: c = p.new_variable(dim=2)

sage: # relaxed value
sage: r = p.new_variable(dim=2)

sage: E = lambda x,y : (x,y) if x<y else (y,x)

sage: MAD = 1-1/(Integer(g.order())*2)

sage: # Partition of the edges
sage: for u,v in g.edges(labels=None):
...     p.add_constraint(sum([c[i][E(u,v)] for i in range(k)]), max=1, min=1)

sage: for i in range(k):
...     # r greater than c
...     for u,v in g.edges(labels=None):
...         p.add_constraint(r[i][(u,v)] + r[i][(v,u)] - c[i][E(u,v)], max=0, min=0)
...         # Maximum degree 2
...         for u in g.vertices():
...             p.add_constraint(sum([c[i][E(u,v)] for v in g.neighbors(u)]),max = 2)
...             # no cycles
...             p.add_constraint(sum([r[i][(u,v)] for v in g.neighbors(u)]),max = MAD)

sage: p.set_objective(None)
sage: p.set_binary(c)

sage: c = p.get_values(c)

sage: gg = g.copy()
sage: gg.delete_edges(g.edges())
sage: answer = [gg.copy() for i in range(k)]
sage: add = lambda (u,v),i : answer[i].add_edge((u,v))

sage: for i in range(k):
...         for u,v in g.edges(labels=None):
...             if c[i][E(u,v)] == 1:
...                 add((u,v),i)
\end{lstlisting}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{H-minor}

For more information on minor theory, please see \\\url{http://en.wikipedia.org/wiki/Minor_\%28graph_theory\%29}\\ It is a wonderful subject, and I do not want to begin talking about it when I know I couldn't freely fill pages with remarks :-)

For our purposes, we will just say that finding a minor $H$ in a graph $G$, consists in :
\begin{enumerate}
\item Associating to each vertex $h\in H$ a set $S_h$ of representants in $H$, different vertices $h$ having disjoints representative sets
\item Ensuring that each of these sets is connected (can be contracted)
\item If there is an edge between $h_1$ and $h_2$ in $H$, there must be an edge between the corresponding representative sets
\end{enumerate}

Here is how we will address these constraints :

\begin{enumerate}
\item Easy
\item For any $h$, we can find a spanning tree in $S_h$ (an acyclic set of $|S_h|-1$ edges)
\item This one is very costly.

To each {\em directed} edge $g_1g_2$ (I consider $g_1g_2$ and $g_2g_1$ as different) and every edge $h_1h_2$ is associated a binary variable which can be equal to one only if $g_1$ represents $h_1$ and $g_2$ represents $g_2$. We then sum all these variables to be sure there is at least one edge from one set to the other.
\end{enumerate}

\newpage
{\bf LP Formulation :}
\begin{itemize}
\item Maximize : nothing
\item Such that :
  \begin{itemize}
  \item A vertex $g\in V(G)$ can represent at most one vertex $h\in V(H)$
    $$\forall g\in V(G), \sum_{h\in V(H)}rs_{h,g}\leq 1$$
  \item An edge $e$ can only belong to the tree of $h$ if both its endpoints represent $h$
    $$\forall e=g_1g_2\in E(G), t_{e,h}\leq rs_{h,g_1}\text{ and }t_{e,h}\leq rs_{h,g_2}$$
  \item In each representative set, the number of vertices is one more than the number of edges in the corresponding tree
    $$\forall h, \sum_{g\in V(G)}rs_{h,g}-\sum_{e\in E(G)}t_{e,h} = 1$$
  \item No cycles in the union of the spanning trees
    \begin{itemize}
    \item Each edge sends a flow of 2 if it is taken
      $$\forall e=uv\in E(G), x_{e,u} + x_{e,v} = 2\sum_{h\in V(H)}t_{e,h}$$
    \item Vertices receive strictly less than 2
      $$ \forall v\in V(G), \sum_{e\in E(G)\atop e\sim v}x_{e,k,v}\leq 2-\frac 2 {|V(G)|}$$
    \end{itemize}
  \item $arc_{(g_1,g_2),(h_1,h_2)}$ can only be equal to 1 if $g_1g_2$ is leaving the representative set of $h_1$ to enter the one of $h_2$. (note that this constraints has to be written both for $g_1, g_2$, and then for $g_2, g_1$)
    $$\forall g_1, g_2\in V(G), g_1\neq g_2, \forall h_1h_2\in E(H)$$$$arc_{(g_1,g_2),(h_1,h_2)}\leq rs_{h_1, g_1} \text{ and }arc_{(g_1,g_2),(h_1,h_2)}\leq rs_{h_2, g_2}$$
  \item We have the necessary edges between the representative sets
    $$\forall h_1h_2\in E(H)$$$$\sum_{\forall g_1, g_2\in V(G), g_1\neq g_2}arc_{(g_1,g_2),(h_1,h_2)}\geq 1$$
  \end{itemize}
\item Variables
  \begin{itemize}
  \item $rs_{h,g}$ binary (does $g$ represent $h$ ? rs = ``representative set'')
  \item $t_{e,h}$ binary (does $e$ belong to the spanning tree of the set representing $h$ ?)
  \item $x_{e,v}$ real positive (flow sent from edge $e$ to vertex $v$)
  \item $arc_{(g_1,g_2),(h_1,h_2)}$ binary (is edge $g_1g_2$ leaving the representative set of $h_1$ to enter the one of $h_2$ ?)
  \end{itemize}
\end{itemize}

Here is the corresponding Sage code:

\begin{lstlisting}
sage: g = graphs.PetersenGraph()
sage: H = graphs.CompleteGraph(4)

sage: p = MixedIntegerLinearProgram()

sage: # sorts an edge
sage: S = lambda (x,y) : (x,y) if x<y else (y,x)

sage: # rs = Representative set of a vertex
sage: # for h in H, v in G is such that rs[h][v] == 1 if and only if v
sage: # is a representant of h in g
sage: rs = p.new_variable(dim=2)

sage: for v in g:
...         p.add_constraint(sum([rs[h][v] for h in H]), max = 1)

sage: # We ensure that the set of representatives of a
sage: # vertex h contains a tree, and thus is connected

sage: # edges represents the edges of the tree
sage: edges = p.new_variable(dim = 2)

sage: # there can be a edge for h between two vertices
sage: # only if those vertices represent h
sage: for u,v in g.edges(labels=None):
...         for h in H:
...             p.add_constraint(edges[h][S((u,v))] - rs[h][u], max = 0 )
...             p.add_constraint(edges[h][S((u,v))] - rs[h][v], max = 0 )

sage: # The number of edges of the tree in h is exactly the cardinal
sage: # of its representative set minus 1

sage: for h in H:
...         p.add_constraint(
...           sum([edges[h][S(e)] for e in g.edges(labels=None)])
...           -sum([rs[h][v] for v in g])
...           ==1 )

sage: # a tree  has no cycle
sage: epsilon = 1/(5*Integer(g.order()))
sage: r_edges = p.new_variable(dim=2)

sage: for h in H:
...         for u,v in g.edges(labels=None):
...             p.add_constraint(
...              r_edges[h][(u,v)] + r_edges[h][(v,u)] >= edges[h][S((u,v))])
...         for v in g:
...             p.add_constraint(
...                sum([r_edges[h][(u,v)] for u in g.neighbors(v)]) <= 1-epsilon)

sage: # Once the representative sets are described, we must ensure
sage: # there are arcs corresponding to those of H between them
sage: h_edges = p.new_variable(dim=2)

sage: for h1, h2 in H.edges(labels=None):
...         for v1, v2 in g.edges(labels=None):
...             p.add_constraint(h_edges[(h1,h2)][S((v1,v2))] - rs[h2][v2], max = 0)
...             p.add_constraint(h_edges[(h1,h2)][S((v1,v2))] - rs[h1][v1], max = 0)
...             p.add_constraint(h_edges[(h2,h1)][S((v1,v2))] - rs[h1][v2], max = 0)
...             p.add_constraint(h_edges[(h2,h1)][S((v1,v2))] - rs[h2][v1], max = 0)



sage: p.set_binary(rs)
sage: p.set_binary(edges)
sage: p.set_objective(None)
sage: p.solve()                                         # optional - GLPK,CBC,CPLEX
0.0

sage: # We can now build the solution found as a
sage: # dictionary associating to each vertex of H
sage: # the corresponding set of vertices in G

sage: rs = p.get_values(rs)

sage: from sage.sets.set import Set
sage: rs_dict = {}
sage: for h in H:
...         rs_dict[h] = [v for v in g if rs[h][v]==1]
\end{lstlisting}
